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Abstract .  A  variant  of  the  classical  pseudo-arclength 

continuation  method  is  proposed.  Basically,  the  method  can 

‘j 

be  viewed  as  pseudo-arclength  continuation  in  (r, A) £ space 
where  r  is  a  functional  of  the  solution.  Another  difference 
is  a  three-parameter  predictor  instead  of  the  standard  Euler 
step.  This  predictor,  as  well  as  the  Newton  corrector 
iteration,  are  justified  and  some  numerical  results  for 
reaction-diffusion  equations  are  presented.  The  method 
provides  a  simple  algebraic  check  for  symmetry-breaking 
bifurcation,  the  most  common  type  of  secondary  bifurcation 
in  physical  examples. 
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Introduction 


Boundary  value  problems  for  partial  differential  equations 
that  describe  phenomena  in  many  fields  of  application  are 
usually  nonlinear  and  depend  on  one  or  several  parameters  of 
the  underlying  system.  These  parameters  enter  in  various  ways. 
The  coefficients  or  right-hand  sides  of  the  differential 
equations  may  depend  on  them,  the  boundary  conditions  or,  for 
example,  the  domain.  With  respect  to  these  parameters,  in  gen¬ 
eral,  bifurcation  and  non-unicity  effects  occur.  For  a  list  of 
different  problems  from  the  applications  as  well  as  a  first 
introduction  to  analytic  and  numerical  aspects  we  refer  to  [8], 
A  survey  of  the  state  of  the  art  may  be  obtained  from  [7], 

Frequently,  the  only  way  to  explore  the  solution  manifolds 
of  these  problems  is  to  continue  along  the  solution  branches, 
in  the  case  of  multi-parameter  dependence,  with  respect  to  one 
of  the  parameters  for  fixed  values  of  the  others.  This 
approach  yields  a  quan t i t a t i ve ly  satisfactory  approximat  ion  if, 
for  example,  suitable  discretizations  are  used,  while  it  may  be 
difficult  to  get  a  full  qualitative  impression  in  this  way. 
Other  analytic  techniques  yielding  the  latter  but  not  also  a 
quantitatively  useful  result,  may  well  supplement  the  numerical 
methods . 

For  the  continuation  along  solution  branches,  different 
methods  have  been  proposed  and  applied.  We  restrict  considera¬ 
tions  here  to  the  methods  introducing  an  arclength- like 
additional  parameter.  One  of  the  earlier  and  frequently 
referenced  papers  is  [6],  To  be  more  specific,  let  the 


parameter-dependent  nonlinear  problem  be  given  by 

(1.1)  G ( u, a)  =  0 

where  u  is  an  element  of  a  Banach  space  X  and  G  maps  XxR^  into 
X.  We  denote  by  X  the  component  of  at  that  is  presently  not 
fixed  and,  for  simplicity,  suppress  the  others.  The  pseudo- 
arclength  method  of  [6]  then  introduces  an  additional  parameter 
a  and  augments  (1.1)  by  a  normalizing  condition 

G (u , X )  =  0, 

(1.2) 

N  (  u(  a) ,  X(  o  )  ,  o)  =  0  . 

Thus,  u  as  well  as  X  are  locally  parametrized  by  o.  In  the 
following  section,  we  will  Introduce  a  slight  modification  of 
this  approach.  In  applications,  one  is  frequently  interested 
in  a  functional  of  u,  for  example,  a  norm,  the  value  of  the 
definite  integral  of  u  over  some  region,  etc.  This  quantity 
may  represent  a  physically  important  parameter  just  as  the 
components  of  a.  For  graphical  representations  also,  usually 
such  a  functional  is  depicted  versus  one  or  two  of  the  a's. 
Using  a  functional  of  u,  therefore,  is  not  equivalent  to 
introducing  another  artificial  parameter. 

After  introducing  the  basic  continuation  method,  we 
discuss  a  new  predictor  to  be  used  in  conjunction  with  it  in 
Section  3.  An  analysis  of  the  corrector  iteration  is  given  in 
the  following  section,  and  in  the  last  section,  some  numerical 
results  are  presented. 
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2 .  The  Continuation  Method 
For  the  problem 

(2.1 )  G(u, X)  =  0,  G : XxR  +  X 

let  s  denote  the  arclength  along  a  solution  arc  (u(s),A(s)). 
Under  sufficient  smoothness  and  regularity  assumptions  this 
latter  dependence  is  differentiable  and  thus  (G°  =  Gu(uo»Aq)) 

(2.2a)  G°  u0  +  G?  X0  =  0  , 

U  A 

•  2  *2 

(2.2b)  llu0»  +  X0  =  1 

in  a  solution  (uq,X0)  =  (  u ( sq ) , A ( so  )  )  of  (2.1).  Here,  sub¬ 
scripts  stand  for  partial  derivatives  and  a  dot  for  differen¬ 
tiation  with  respect  to  s. 

In  [6]  it  was  proposed  to  augment  (2.1)  by  the  following 
normalizing  condition: 

(2.3)  No(u,X,o)  =  9  uo(u-uo)  +  ( 2  -  0  )  X 0 ( A -A 0  )  -  6o  =  0 

which  approximates  (2.2b).  o  is  thus  called  the  pseudo- 
arclength  parameter.  6c  is  the  steplength  in  c.  0  satisfying 
0  <  0  <  2  was  meant  to  be  a  weighting  parameter. 

Let  r(u)  denote  a  functional  of  u.  For  simplicity,  let 

(2.4)  r  (  u  )  =  Hull 


where  u • J  is  the  norm  of  u  which  we  assume  to  be  differenti¬ 
able.  Other  functionals  may  be  chosen,  as  for  example,  in  VLSI 
problems  ([2]),  but  in  the  following  we  restrict  considerat ions 
to  (2.4).  Quite  frequently,  this  parameter  is  used  in  a  graph¬ 
ical  representation  of  the  solution  manifold  and  also  has  a 
physical  s ign i f icance . 

The  following  augmenting  equation  leads  to  a  pseudo- 
arclength  method  in  (r,X)-space. 

(2.5)  N  r  (  u ,  X  ,  o )  5  9  r  g  (  r  -  r  g  )  +  (2-8)X0(X-X0)  -  6a  =  0, 

where  0  <  8  <  2,  rg  =  3ugll,  rg  =  (—  Hull) 

da  0  =  °0 

While  the  classical  pseudo-arclength  method  with  Euler 

predictor  using  (2.3)  is  in  general  optimal  in  a  local  sense 

among  continuation  procedures  utilizing  only  first  derivatives 

and  depending  only  on  one  parameter,  in  this  case  the  step- 

length  along  the  tangent,  (2.5)  will  be  combined  with  a  more 

exp'nsive  predictor  to  yield  in  many  examples  a  better  global 

continuation  behaviour.  This  will  be  demonstrated  in  5. 

Let  us  assume  from  now  on  that  X  is  a  Hilbert  space, 

1  /  2 

sufficient  regularity  is  given,  and  H • J  =  (•,•)  ,  then 

(2.6)  r g  =  (ug, Ug)  /rg  . 

From  (2.5)  we  see  that  the  three  values  0,1,2  for  9  are  special 
in  the  sense  that  for  9=0  (2)  the  augmenting  equation 
characterizes  points  with  a  fixed  X(r)-value,  while  for  9  =  1 


the  points  lie  on  a  hyperplane  orthogonal  to  the  solution  arc 
in  ( r  o,  A0). 

This  suggests  to  continue  to  target  values  in  the  two 
physical  parameters  r  and  A.  Such  a  strategy  seems  appropriate 
as  long  as  starting  from  a  point  on  the  solution  curve  another 
point  corresponding  to  a  suitably  picked  target  value  can  be 
computed  in  a  few  number  of  (corrector)  iterations.  This  will, 
of  course,  depend  on  the  starting  guess  and  thus  on  the  method 
used  for  predicting  this  point.  Only  the  combination  of  (2.5) 
with  the  predictor  presented  in  the  next  section  makes  this 
continuation  the  effective  and  powerful  algorithm  it  has  proven 
to  be  at  least  for  a  series  of  problems  from  the  applications. 


3 .  The  Predictor 

In  continuing  from  a  point  (uo>Ao)  on  a  solution  branch 
most  continuation  methods  first  generate  a  starting  guess 


and  then  iteratively  determine  a  new  solution  point. 


The  pseudo-arclength  method  in  its  standard  version  uses  a 


first  order  Euler  predictor 


(uo  + 

A  o  + 

Subsequently,  the  augmented  system  (2.1),  (2.3)  is  iteratively 

solved  ('corrector').  The  computed  solution  for  0  =  1  is  thus 
located  in  the  intersection  of  the  solution  manifold  and  a 

•  •  J 

hyperplane  perpendicular  to  the  tangent  vector  (  u  q  ,  A  0  )  at  the 
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previous  point  and  a  distance  <5o  from  this  point.  It  is  not 
immediately  apparent  how  the  parameter  9  would  have  to  be 
chosen  to  improve  the  quality  of  the  predicted  point. 

It  is  clear  that  with  respect  to  the  information  used, 

•  • 

namely  u0  and  X0,  the  Euler  predictor  step  is  locally,  i.e., 
for  6o  0,  optimal,  but  that,  for  example,  in  highly  curved 
parts  of  the  solution  arc,  only  a  rather  small  So  will  lead  to 
convergence  of  the  corrector  iteration. 

In  order  to  improve  the  quality  of  the  predictor  either  a 
higher  order  method  making  use  of  higher  derivatives  of  G  may 
be  used  or  in  some  other  way,  it  has  to  be  assured  that  the 
predicted  point  better  approximates  a  solution  of  the  nonlinear 
eigenvalue  problem.  The  first  approach  is  suitable  in  particu¬ 
lar  if  the  nonlinearities  have  a  simple  form  as,  for  example, 
low  order  polynomials.  Here,  we  propose  to  make  use  of  only 
the  first  order  information  also  used  in  the  Euler  predictor. 
Instead  of  a  single  parameter,  three  are  used  in 


In  case  a  *  S,  y  =  0  the  predicted  point  is  in  a  no n - 1 ange n t  ia  1 
direction  from  (uo,Xq).  This  is  slightly  more  general  than  (3.1) 
but,  of  course,  treats  X  as  a  special  variable  as  opposed  to 
methods  that,  for  a  f i n i t e - d imens i on  a  1  problem,  give  X  no  prefer¬ 
ential  treatment  (  [ 1 2  ] ) .  The  justification  for  (3.2)  can  only  be 
that  X,  in  fact,  represents  an  important  physical  parameter. 


,m\  A  wflfc  mS  m, 


More  explanation  is  needed  for  the  coefficient  of  uo  in 

(3.2) .  Before  this  will  be  done,  we  complete  the  definition  of 
the  predictor.  A  system  of  three  equations  is  used  to 
determine  the  parameters  a,  3  and  y. 

N(up,Xp,o)  =  0  , 

(3.3)  ( u0 ,G( Up , A  ) )  =  0  , 

( u0 ,  G( Up  ,  Ap  )  )  =  0  . 

In  addition  to  the  normalizing  equation,  two  orthogonality 
relations  with  respect  to  the  scalar  product  on  X  are  imposed. 
The  basis  for  these  equations  is  twofold.  First,  it  assures 
that  the  predicted  point  is  a  weak  solution  of  (2.1)  with 
respect  to  the  subspace  spanned  by  uo  and  uo .  On  the  other 
hand,  (3.2)  shows  that 

(3.4)  (u  ,  G(u  ,A  ))  =  0 

P  P  P 

holds,  which  defines  a  generalized  Ray  1 e i gh - quo t  ien t . 

The  system  (3.3)  in  general  need  not  have  a  solution 
(a,3,y)  or  such  a  solution  is  not  necessarily  unique.  If,  for 
example,  ug  and  ug  are  parallel,  then  the  rank  of  the  Jacobian 
of  (3.3)  is  at  most  two.  Therefore,  a  least  squares  solution 
of  (3.3)  is  to  be  computed.  In  order  to  be  able  to  analyze 
several  important  features  of  the  continuation  method,  we  have 
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to  more  closely  consider  the  case  of  bifurcation  points.  Let 
us  introduce  the  notation 


(G(u,A)  \ 

)  ,  y  =  (u,x) 

N  (  u  ,  A  ,  o  )  J 

and  use  subscripts  for  partial  derivatives.  (ug,A0)  is  a 
simple  pitchfork  (symmetric)  bifurcation  point  of  (2.1)  if  the 
following  conditions  hold 

N(G°)  =  span  { 4> 0 }  ,  Ifo11  =  1  , 

R(G°)  is  closed  and  codim  R(G°)  =  1  , 

u  u 

(3.6) 

4>oG°  =  0  ,  *}>  oG  uU^0  4>0  =  0  > 

b  -  +  GuUvO<i>o)  *  0 

where  vq  is  the  unique  solution  of 
(3.7)  C°v0  +  cj  =  0,  4>  o v  o  =  0  . 

N,  R  denote  nullspace  and  range,  respectively,  and 
N(G°  )  =  span  J4>0}  ,  4>  0  4>  o  =  1  • 


From  differentiating  (2.1)  twice  with  respect  to  s  it  follows 
that  (cf.  [**•])  the  algebraic  bifurcation  equation 


(3.8) 


(cXq  +  2  b  ?  i  )  X  q  =  0 


holds  where  c  =  <J>  o  ( G 0  vqvo  +  2G°,vq  +  G?,)  and  C2  +  2X2  =  1. 

UU  UA  AA  j  g 

The  two  solution  branches  passing  through  (uo,Xq)  are  thus 
given  by 


(3.9) 


•  n 

u(s)  =  uq  +  (  s-so  )  (  X3  vq  +  ^  1 4>  0  )  +  0((s-so)  )  , 
X  ( s )  =  X  q  +  ( s - s  o  )  Xq  +  0((s-sq)2) 


These  pitchfork  bifurcation  points  are  not  generic  in  the  sense 
that  they  are  in  general  not  present  under  perturbations  of 
the  problem  as,  for  example,  introduced  by  a  discretization. 

An  even  further  specialized  class  of  bifurcation  points,  the 
symmetry-breaking  bifurcation  points,  however,  have  this 
property . 

Let  us  assume  for  simplicity  that  S  is  a  linear  operator 
on  X,  S2  =  I,  S  *  I  such  that 

(3.10)  G(Su,X)  =  SG(u,X),  for  all  u  e  X,  X  e  R  . 


Then  X  may  be  split  as 
(3.11  )  X  =  Xs©  Xa  , 

where  X  ^  =  {u  e  X,  u  =  Su }  contains  the  symmetric  elements  with 
respect  to  S  and  X  =  ju  e  X,  u  =  -Su}  ,  the  anti-symmetric 
elements.  The  pitchfork  bifurcation  point  (ug,X0)  is  then 
sy mme t ry - br eak i ng  if 


(3.12)  uq  £  X  ,  v0  e  X  ,  ♦o  e  X 

S  S  a 


On  the  branch  corresponding  to  the  solution  Xg  =  0  of  (3.8)  it 
thus  holds  that 


(3.13)  u0  =  —  tO  -  X 


while  on  the  other  branch  c  =  0  implies  ug  e  .  On  the 
symmet ry- break ing  branch  we  have  from  (2.6)  if  S  =  S* 


(3.14)  X0  =  0  ,  r o  =  0 


and  thus,  in  this  case,  Nf  in  (2.5)  is  not  well  defined.  The 
restriction  to  the  two-dimensional  ( r , X  ) -subspace  introduces 
additional  singularities  that  are  not  present  in  the  classical 
pseudo-arclength  method.  A  natural  modification  is  to  intro¬ 
duce  as  a  parameter  the  norm  of  the  anti-symmetric  part  u  of  u 
corresponding  to  the  splitting  (3.11).  This  was  done  explicit¬ 
ly  in  [lO].  Since  on  the  branch  of  nonsymmetric  solutions 
u  =  uq  e  X  ,  asymptot ical ly  in  the  bifurcation  point  from 
(3.9),  we  propose  instead  in  the  general  case  to  use  N  =  in 
(3.3)  whenever  J  r  q  |  /  il  u  q  tl  is  small. 

Let  us  assume  now  that  starting  at  (ug,X0)  a  point  (u,X) 
on  the  s  ymme  t  r  y -b  reak  i  ng  branch  with  r  =  Hull  <  rg  =  »  ug  i 
is  to  be  computed.  For  y  =  0  (3.2),  (3.12),  (3.13)  imply  that 


always  :l >  rg 


This  is  one  of  the  reasons  that  y  was 
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introduced  in  (3.2).  Another  reason  is  that  in  the  linear 
eigenvalue  problem,  computations  along  the  eigenfunction 
branches  require  a  renormalization  only  and  are  not  prone  to 


cancellation  if  in  the  predictor  step  u^  is  actually  obtained 


by  renormalization,  i.e.,  6  =  0  in  (3.2). 

Let  us  summarize  the  normalizing  condition  that  is 
proposed  in  the  various  cases. 


j  r  0  j  <  e  II  u  o  »  : 
(3.15)  |  r  o  |  >  ( 1  -  e  ) «  u  o  ll  : 


N  =  N  ,  9  *  0 
a  ’ 


N  =  Nr,9  *  0,  B  =  0  , 


e  II  u  o  ll  £  r0  £  ( 1  - e  )  II  u j,  II  :  N  =  Nr,  0  £  9  £  2 


In  the  following  section,  we  will  address  the  question  of 
the  regularity  of  the  corrector  Dacobian  in  the  various  cases 
of  (3.15). 


4.  The  Corrector 


The  point  (up>*p)  computed  by  the  predictor  described  in 


the  previous  section  is  used  as  usual  as  starting  guess  for  an 
iterative  solution  of  the  augmented  system  (3.5).  This  is  done 
by  Newton's  method  and  thus  the  Dacobian 


(4.1  ) 
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has  to  be  evaluated  and  inverted.  The  following  lemma  is  a 
simplified  version  of  Lemma  2.8  in  [6]. 


Lemma  4 .  1  Let  X  be  a  Banach  space  and  consider  the  linear 
operator  As  XxR  XxR  of  th$  for?!! 


where  A :  X  X ,  B :  R  +  X ,  C* :  X  ♦  R ,  L) :  R  +  R  . 

(i)  If  A  is  nonsingular  then  A  is  nonsingular  iff 

(4.2)  E  =  D  -  C*A_1B  is  nonzero. 

(ii)  If  A  is  singular  and  dim  N  ( A )  =  codim  R(A)  =  1  then 
A  is  nonsingular  iff 


(4.3) 


B  \  R ( A  )  ,  C*  f  R ( A  *  )  . 


Proposition  4.2  For  an  appropriate  choice  of  the  continuation 
parameter,  i.e.  X  respectively  r,  in  (4.1)  is  regular 
except  at  bifurcation  points. 


Proof  The  quantity  E  In  (4.2)  is  easily  computed  as 

*  •  2  • 

E  =  ( 2 -  0 ) X  g  +  Gto/Xq  if  n  =  N  .  If  we  accept  as  an  appropri¬ 
ate  choice  of  the  continuation  parameter  not  to  pick  r  respec¬ 
tively  X  to  continue  past  a  turning  point  in  r  respectively  X 
then  we  have  0  *  0  if  X0  =  0  respectively  0  *  2  if  rg  =  0,  so 


that  Fy  is  regular  in  regular  points  for  N  =  N^.  In  bifurca¬ 
tion  points  e  FKG^)  so  that  is  always  singular  from  (4.3) 
for  N  =  Nr  and  N  =  N  .  It  remains  to  consider  turning  points. 

If  X0  =  0  and  G,  i  R(G  )  then  the  last  condition  in  (4.3) 

A  U 

.  .  -2 

is  equivalent  to  (  erouo/ro  ,  uo  )  =  tiro  *  0  .  So  i  f  6  *■  0 ,  rg  *  0 

then  Fy  is  regular  for  N  =  unless  ro  =  0.  This  situation 

occurs  in  symmetry-breaking  pitchfork  bifurcation  points. 

These  points  appear  to  be  the  most  common  type  of  secondary 

bifurcation  in  physical  examples.  Finite  dimensional  problems 

exhibiting  secondary  bifurcation  often  arise  from  symmetry 

preserving  di scr e t iza t i ons  of  such  continuous  problems.  It  is, 

therefore,  important  that  the  continuation  method  can  handle 

these  points.  In  addition  to  that,  the  proposed  algorithm 

allows  to  detect  these  points  by  checking  (3.14).  We  refer  to 

Example  1  in  [  1 1  ]  where  two  symmetry-breaking  bifurcation 

points  are  computed  and  to  the  second  example  in  5.  According 

to  (3.13)  N  =  N  is  chosen  in  these  points  since  N  is  not 
o  r  r 

well  defined.  If  the  point  is  not  a  bifurcation  point  both 
conditions  of  (4.3)  are  satisfied  and  F^  is  again  regular. 

The  above  proposition  does  not  specify  in  any  more  detail 
what  happens  in  bifurcation  points.  If  we  consider  F"1  in  the 
neighbourhood  of  a  simple  bifurcation  point  v  =  y(so)  and 
measure  the  growth  of  JJ F ~  1  )i  as  a  function  of  (s-s0)_1  as  in  [4] 
then  we  obtain  the  following  result  along  the  lines  of  proofs 
in  [4,3]. 


■fi 
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Proposition  4 . 3  If  the  normalization  condition  M  is  chosen 
according  to  (3.15)  then  on  any  of  the  branches  passing  through 
a  simple  bifurcation  point  yg  =  ( u  ( so  ) »  Msg))  of  (2.1)  for 
6  >  0  sufficiently  small 

(4.4)  » F  ~  1  ( y  (  s  ) )  II  =  0(  |  s  -  s  o  j  ~  1  )  ,  0  <  js-s0j  <  6  . 

The  growth  behaviour  of  (4.4)  allows  the  application  of  the 
Newton-Kantorovich  theorem  to  show  convergence  for  the  starting 
point  obtained  by  the  Euler  predictor  ([4]).  For  the  combina¬ 
tion  with  the  predictor  (3.2),  (3.3)  a  different  analysis  would 

have  to  be  done.  One  of  the  facts  that  complicate  such  an 
analysis  considerably  is  the  least  squares  solution  of  (3.3). 

We  do  not  present  such  a  theory  here  and  refer  to  an  illustra¬ 
tion  of  the  practical  merits  of  the  method  by  the  numerical 
results  in  [ 1 1 ]  and  in  the  following  section. 

5 .  Numerical  Results 

We  present  here  only  a  demonstration  of  the  behaviour  of 
the  proposed  continuation  method.  The  method  was  developed 
and  implemented  in  the  program  PLTKG  [l]  in  cooperation  with 
R .  Bank.  The  class  of  second  order  nonlinear  parameter- 
dependent  boundary  value  problems  solved  by  PLTMG  includes  many 
interesting  app  1  i ca t  ions  . 

It  is  of  interest  how  the  continuation  method  proposed 
above  compares  to  standard  methods  as,  for  example,  the 


classical  pseudo-arclength  method,  8  =  1  in  (2.3).  The  gener¬ 
alized  inverse  iteration  ([9])  that  strongly  influenced  the 
method  presented  here  was  compared  to  different  continuation 
methods  in  [9]  and  to  a  predecessor  of  the  present  method  using 
a  two-parameter  predictor  and  (2.3)  in  [ 1 0 ]  .  The  latter  paper 
also  contains  comparisons  with  a  variant  of  the  generalized 
inverse  iteration  in  the  neighbourhood  of  symmetry-breaking 
bifurcation  points. 

An  extensive  comparison  of  the  present  method  with  the 
pseudo-arclength  method  is  beyond  the  scope  of  this  paper.  We 
restrict  the  consideration  to  one  example  exhibiting  turning 
points  and  thereby  illustrate  the  advantages  of  a  multi¬ 
parameter  predictor.  The  second  example  has  a  more  complicated 
solution  structure  and  serves  to  demonstrate  the  computation  of 
singular  points,  the  switching  of  branches  at  bifurcation 
points,  the  detection  of  s ymme t ry -b reak i ng  bifurcation  points 
and  the  switching  of  continuation  parameters  in  a  multi- 
parameter  problem. 

The  first  example  is  the  two-parameter  Bratu  problem 
(5.1)  Au  +  Aexp(u/(1+su))  in  ft 

on  the  unit  square  ft  with  homogeneous  Dirichlet  boundary  condi¬ 
tions.  This  problem  was  solved  frequently  (see,  for  example, 
[3,1l]).  For  e  =  .1  the  solution  curve  has  two  turning  points 
and  we  present  results  for  continuing  up  to  and  beyond  the 


L«, 


first  one.  We  point  out  that  this  turning  point  is  far  from 
being  nearly  angular  and  therefore  does  not  cause  particular 
difficulties  for  the  pseudo-arclength  method.  For  a  case  of  a 
nearly  angular  turning  point,  see  Example  2  in  [ll]. 

On  a  standard  triangulation  of  ft,  with  25  vertices,  we 
have  discretized  (5.1)  by  linear  finite  elements.  The  result¬ 
ing  f in i t e - d ime ns i ona 1  system  was  first  solved  by  the  pseudo- 
arclength  method  starting  at  the  origin.  We  assumed  that  the 
step-picking  procedure  does  not  cut  back  the  steplength  a  in 
the  Euler  predictor,  i.e.,  a  =  B,y  =  0  in  (3.2),  enough  near 
the  turning  point  to  avoid  failure  of  the  corrector  iteration. 
In  this  latter  case,  the  steplength  is  cut  back  by  .5  until  the 
corrector  converges.  In  the  last  column  of  Table  5.1  the  num¬ 
ber  of  corrector  iterations  is  listed  and  an  asterisk  denotes 
failure  to  converge  after  the  given  number  of  iterations  while 
two  asterisks  denote  failure  to  yield  a  descent  direction  for 
the  damped  Newton  process. 

Subsequently,  we  have  used  our  continuation  method  to 
continue  from  each  of  the  points  in  Table  5.1  below  the  turning 
point  directly  to  the  point  above  that  with  " u "  =  2.61  which 
was  finally  reached.  The  results  in  Table  5.2  show  that  at 
most  3  corrector  steps  are  needed.  The  three  predictor  parame¬ 
ters  are  also  given  to  show  for  this  example  the  effect  of 
introducing  y  and  the  early  'detection'  of  the  turning  point 


Table  5.2  Continuation  with  proposed  method  for  each  point 
below  turning  point  in  Table  5.1  directly  to 
solution  with  norm  r  =  2.61 

A  suitable  step-picking  procedure  may  be  able  to  reduce  the 
number  of  failing  corrector  Iterations  for  the  pseudo- arc leng t h 
method.  But  the  fact  remains  that  several  small  steps  are 
necessary  to  overcome  the  turning  point.  The  next  example  has 
a  more  complex  solution  diagram  including  bifurcation.  It 
illustrates  several  other  features  of  the  continuation  method 
and  its  implementation. 

Let  1 f  be  a  domain  that  is  obtained  by  connecting  two 
regions  in  the  left  and  right  half-plane  situated  symmetrically 
w.r.t.  the  y-axis  by  a  channel  of  width  e  such  that  Q  is 
symmetric,  too.  Homogeneous  Neumann  boundary  conditions  are 
prescribed  along 


3u 


2  3 


I 


au  -  u 


(5.1) 


A  u  +  Xu 


in  n 


The  applications  of  this  react  ion- di f f us  ion  equation  include  in 
particular  a  selection-migration  model  from  population  genetics 
(cf.[l3]).  We  are  interested  in  the  stationary  solutions  only. 
The  bifurcation  behaviour  of  (5.1)  w.r.t.  the  three  parameters 
e,A,a  is  given  in  [  1 3  ]  .  For  a  =  0,  for  example,  secondary 
bifurcation  takes  place  from  the  branch  bifurcating  from  the 
trivial  solution  at  the  first  nonzero  eigenvalue  A£  of  the 
linearization  of  (5.1).  For  a  *  0,  and  for  this  case  some 
results  will  be  presented  below,  secondary  symmetry-breaking 
bifurcation  occurs  from  the  branch  of  constant  solutions 
bifurcating  at  the  origin. 

Because  of  symmetry  and  the  homogeneous  Neumann  boundary 
condition,  it  suffices  to  take  as  a  non-convex,  say, 

L-shaped  domain  with  vertices  (0,0),  (1,0),  ( 1 , . 5  )  ,  (.5, .5), 
(.5,1),  (0,1).  Thus,  e  is  fixed  and  the  parameters  a  and  A  are 
left.  For  a  >  0  the  following  schematic  bifurcation  diagrams 
were  obtained  in  [ 1 3 ]  for  a  point  value  of  u  versus  A. 


The  elliptic  version  of  (5.1)  was  discretized  on  a 
relatively  coarse  triangulation  with  21  vertices  and  linear 
finite  elements  were  used  to  compute  the  approximate  solution 
u^  .  Better  approximations  at  specified  points  of  the  solution 
branches  may  be  obtained  through  PLTMG  by  a  multi-grid  method. 

The  first  nontrivial  eigenvalue  was  X£h  =  6.73.  In  fact, 
for  a  =  5  no  secondary  bifurcation  was  observed.  For  a  =  S  the 
two  secondary  bifurcations  points  in  Figure  5.1b  are  located  at 
( il,r  i)  =  (-8.17,1.04)  and  (A2»r2)  =  (-14.6,2.42),  while  the 
turning  point  is  at  (*3^3)  =  (-16,3.46). 

We  present  below  the  output  of  PLTMG  for  the  following 
continuation.  Starting  on  the  trivial  solution  in  X  =  0  we 
switch  branches  three  times  using  Method  I  of  [6],  On  the 
branch  of  constant  negative  solutions  we  continue  to  (ij ,ri  ). 
There  we  switch  branches  once  and  continue  on  this  branch  of 
non-constant  solutions  to  (X,r)  =  (-11.5,1.6).  Now  we  keep 
this  X  fixed  and  switch  to  a  as  continuation  parameter.  From 
the  value  a  =  8  we  continue  to  smaller  values.  We  expect  to 
first  come  to  the  constant  solutions  and  continuing  further  to 
reach  a  point  where  no  further  continuation  is  possible.  The 
first  happens  between  a  =  7.5  and  a  =  7.4,  as  a  more  detailed 
computation  shows,  while  the  branch  ends  at  a  =  6.78.  If 
instead  we  continue  in  r,  then  (a,r)  =  (6.78,2.94)  is  a  turning 
point  w.r.t.  a  as  expected. 
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Continue  to  r  = 

1 ,  then  to 

r  =  1.5  and 

check  sign  of 

determinant . 
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0 . 1 00E+0 1 
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Determinant  changes  sign; 

find  bifurcation  point  by 

secant  method 

on  5  . 
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Switch  branches 

0  -0 . 81 7E+01 

and  continue  to  r  =  1.2 

0.104E+01  0.484E-01 

,  r  =  1.6. 
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Switch  continuation  parameter  to  a.  Continue  to  a  = 

7.5,  then  to 

r  =  3  and 

check  sign  of  determinant. 
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0.277E+01 
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-0.989E+00 

-0.21SE+04 
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1  0.678E+01 

0 . 294E  +  0 1 

-0.307E-03 

-0.100E+01 

0.836E+01 

0. 109E-03 

1 

1  0.678E+01 

0.294E+01 

0.844E-05 
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-0 . 232E+00 
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Table  5.1  Continuation  results  for  problem  (5.1) 

In  Table  5.1,  it  denotes  the  number  of  (damped)  Newton  corrector 
iterations,  while  6  is  a  quantity  that  changes  sign  at  simple  bifurcation 

points.  We  point  out  that  after  switching  onto  the  branch  of  non-constant 

•  • 

solutions,  A  and  r  are  both  small,  as  expected,  indicating  a  symmetry¬ 
breaking  bifurcation  point  (cf.(3.14)). 
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